Transition matrix Monte Carlo method for quantum systems 
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We propose an efficient method for Monte Carlo simulation of quantum lattice models. Unlike 
most other quantum Monte Carlo methods, a single run of the proposed method yields the free 
energy and the entropy with high precision for the whole range of temperature. The method is 
based on several recent findings in Monte Carlo techniques, such as the loop algorithm and the 
transition matrix Monte Carlo method. In particular, we derive an exact relation between the DOS 
and the expectation value of the transition probability for quantum systems, which turns out to be 
useful in reducing the statistical errors in various estimates. 

PACS numbers: 02.70.Ss,05.10.Ln,05.30.-d,75.40.Mg 



INTRODUCTION 

The Monte Carlo method for classical and quantum 
lattice models has been improved dramatically since the 
proposal of the Metropolis methodjii]. The introduction 
of the extended ensemble, being one of the ideas that en- 
hanced the Monte Carlo method, is characterized by a 
random walker traveling in the energy space. The most 
vifcU-known methods among the ones based on the ex- 
tended ensemble is the multicanonical method 
this method, one can obtain the density of states (DOS) 
from the histogram of visiting frequency at each value of 
the energy. Since the random walk in the extended en- 
semble methods is typically biased in a complicated way, 
the scaling property of the the traveled distance of the 
walker deviates from that of the free random walk, i.e., 
Rr^Vi- This means that the random- walk nature of the 
method does not necessarily guarantee a quick equilibra- 
tion. However, in many important application, it turns 
out that the random walker can still travel the whole 
range of the energy within a reasonable computational 
time, i.e., a computational time bounded by some poly- 
nomial in the system size. Therefore, the method is quite 
useful for models that has a first-order phase transition 
and/or ones that have local minima in the free energy. 
In addition, the direct calculation of the DOS makes it 
possible to calculate the absolute value of the free energy 
of a given system with high precision. Another advan- 
tage of this type of approach is that thermal averages of 
various quantities can be obtained for the whole range 
of the temperature from only a single run of simulation. 
While the early versions of extended ensemble methods 
still suffered from a slow diffusion, even this problem was 
removed recently by the acceleration method proposed 
by Wang and Landau (WL) 0|, in which the random 
walker is forced to move away whenever it lingers about 
a place for a long time. 

For many important models, classical or quantum, it 



was found that the loop/cluster algorithm is so effi- 
cient that the simulation can be done with no (or negligi- 
ble) difficulty of slowing down at any temperature. How- 
ever, since the algorithm is originally designed for a sim- 
ulation at a fixed temperature, one cannot obtain from a 
single run the whole temperature dependence of various 
quantities nor can one estimate the free energy. There- 
fore, it is natural to wonder if the loop/cluster algorithm 
can be used together with the extended ensemble meth- 
ods. An apparent obstacle is that the energy in the quan- 
tum case is not defined so that constructing a histogram 
on it is straightforward and meaningful. Janke and Kap- 
pler proposed^ that in the cluster algorithm for classical 
spin systems one can use the number of graph elements 
(such as bonds in the Swendsen-Wang algorithm) instead 
of the energy for constructing the histogram, showing a 
way for combining the extended ensemble method and 
the cluster algorithm for classical spin systems. Further- 
more, two of the present authors demonstrated^ the en- 
hancement of the efficiency achieved with the use of the 
WL method, and also suggested that the resulting frame- 
work can be generalized to quantum systems. Later, the 
application to the quantum systems was done with a dif- 
ferent formulation by Troyer, Wessel and Alet |^, who 
also showed the efficiency of the resulting algorithm. 

In the present paper, we propose that a further im- 
provement can be achieved for quantum models by us- 
ing the broad-histogram (BH) relation for estimating the 
DOS. The idea has been partially described in our pre- 
vious paper jl^] for the classical models. For the Ising 
models, it was pointed out that a ratio of the densi- 
ties of states at two neighboring values of the energy can 
be expressed in terms of the microcanonical averages of 
the numbers of distinct potential moves. In the applica- 
tion to the Ising model, the number of potential moves is 
nothing but the number of spins such that flipping one of 
them causes a change in the energy by a given amount. 
Since this is a macroscopic quantity, whereas the number 
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of the visits at each value of the energy is not, we can es- 
timate the DOS more precisely by using the relation than 
we do so directly from the histogram. In the case that we 
discuss here, our DOS is not simply the total number of 
states. Therefore we have to use the generalized form of 
the transition matrix Monte Carlo method • In what 
follows, we describe the method specialized for quantum 
systems. 



THE OUTLINE OF THE ALGORITHM 

The simulation consists of two stages. In the first stage 
we perform short runs using the WL method for obtain- 
ing a working estimate of the DOS. As explained below, 
each run in the first stage is characterized by a parameter 
called the modification factor, which controls the adop- 
tive change of the fictitious weight of the random walk. 
A stronger control is imposed on earlier runs whereas 
a weaker control for the latter runs. Because of this 
adoptively changing weight, the microcanonical ensem- 
ble obtained in each run of the first stage deviates from 
the correct one. This leads to some systematic error. 
However , w e need unbiased estimates of microcanonical 
averages 13] because we use them for obtaining the DOS 
from the BH relation. Therefore, in the second stage we 
perform a single long run with the weight fixed to the one 
obtained in the first stage. This is in contrast to the or- 
dinary WL method where the outcome of the first stage 
is taken as the final estimate of the DOS. The organiza- 
tion of the whole simulation is schematically depicted in 
Fig. ^ We will show this two-stage procedure gives more 
precise results for quantum systems than the one-stage 
procedure of the WL method described in jQi] , within the 
same amount of the total CPU time. 

A little more specifically, the whole simulation can be 
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FIG. 1: The organization of a set of simulation. 



described as follows. We construct a histogram of the 
expansion order or the number of active vertices, which 
we denote as n, instead of the histogram of the energy. 
Accordingly, the fictitious weight f{n) is introduced as a 
function of n, rather than the energy. The Monte Carlo 
dynamics is defined so that the detailed balance holds 
for this weight. We start the first run of the first stage 
by setting the weight f{n) to be constant, i.e., f{n) = 1 
for all n. At the beginning of each run, the counter (or 
the histogram), H{n), of visits to the value n is set to 
be while f{n) is not reset except at the beginning of 
the first run. For the i-tli run we use the modification 
factor Xi = e^^ . Every time the random walker visit a 
state with n{S) ~ n, f{n) is replaced by \if(n) in order 
to penalize the random walker staying at this value of n. 
At the same time, H{n) is increased by one. The i-th 
run is terminated when the histogram becomes flat. To 
be specific, the flatness condition is 

H{n) > 0.8 X ^ H{n ) / ^ 1 



for all n. Since the first stage is only for obtaining a work- 
ing estimate of the DOS, we stop when — InA^ < I0~^ 
whereas in the original WL method — InAi < I0~^ was 
used instead. At the end of the first stage, {f{n))^^ is 
approximately proportional to the generalized DOS dis- 
cussed below. In the second stage, we do the same as 
the first stage but with a longer run and with /(n) fixed 
to be the result of the first stage. In this way, we can 
compute non-biased microcanonical averages of various 
quantities in the second stage. 

In both the first and the second stages, the update of 
the configuration is done using the loop update if the 
update does not involve a change in n. Other similar al- 
gorithms such as the directed loop algorithm can be 
used as well, and the generalization is straight-forward. 
The choice, however, is only of the secondary importance 
for the applications to the models without the magnetic 
field as considered in the present paper. For the loop al- 
gorithm, there are two formulations; one is based on the 
continuous imaginary time and the other is the discrete 
imaginary time. Although the difference in the efficien- 
cies of two formulations is minor, in the present paper we 
have used the discrete time formulation following Troyer, 
Wessel and Alet 9|. In what follows, we reformulate their 
algorithm, use it for deriving the BH relation, and de- 
scribe how the relation is used in the second stage of the 
simulation. 

The systematic error due to the discretization in the 
discrete time formulation can be removed bv starting 
from the high-temperature series expansion rather 
than the path-integral representation. Therefore, we 
start with expressing the partition function Z as a high- 
temperature series expansion truncated at the L-th or- 
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= J2 13- coin), (1) 

Tl = 

where f3 is inverse temperature l/fcsT, /cg is Boltzmann's 
constant, and is a variable that takes on or 1. The 
high-temperature expansion coefficient, is here re- 

garded as the generahzed DOS for the series order n. 
When the Hamiltonian is expressed as a sum of M local 
interaction terms as 7i = X^bli "^b, the generalized DOS 
is given by 

{•n.bi {4>i} l,b 

(2) 

We can simplify the notation by defining a vertex, w, as a 
pair of indices I and b, i.e., v = (l,b). Several other sym- 
bols are also introduced for simplification of the notation: 
r = {i;|7„ = 1}, 4- = (V'i,V'2,---,V'l), "^v = {ipi+ui'i). 

= (V^j+i|(-i76)|V'i), S = (*,r), and n{S) = |r|. 
As a result, we can rewrite the generalized DOS as 

u{n)^ ^ U{S) (3) 

s 

(i(S) = >.) 

where 

^(5) ^ i^i^^ n n ^(*^')- (4) 

Here, A(5'^) = S^i^^^^pi- By comparing these expressions 
to those in Q it is natural to identify n as the number of 
graph elements. Here we call a vertex v for which 7^, = 1 
an active vertex. 

We then replace the factor /3" in ^ by an adjustable 
function /(n), which makes the simulated partition func- 
tion to be 

L 

n=0 S 

with the weight of the state S being 

W{S) = f{n{S))U{S). (5) 

In the first stage, the function fin) is adaptively mod- 
ified so that the probability of having n{S) to be inde- 
pendent of n. In other words, in the first stage we try to 
make f{n)uj{n) constant by adjusting f{n). This goal is 
achieved by the use of the WL method. As a result, we 
obtain u;(n) as the inverse of f{n). 



THE FIRST PHASE AND THE BH RELATION 

The algorithm used for updating S can be viewed as 
the dual Monte Carlo [la, 113 ■ Namely, by regarding the 
variable 71, as a graph element, analogous to the bond 
in the Swendsen-Wang algorithm, one can establish a 
correspondence between the present algorithm and the 
Swendsen-Wang algorithm. One step in a dual Monte 
Carlo simulation consists of two phases; the "graph" F 
is updated in the first phase with "spin configuration" 
^ being fixed and the \1/ is updated in the second phase 
with fixed V. While the first phase changes n, the lat- 
ter does not. In what follows, we describe updating of 
r and show how to estimate w(n) based on estimates of 
macroscopic quantities. The second phase, i.e., how the 
update ^E" with fixed n, is described in the next section. 

In order to relate u){n) to the microcanonical averages 
of some macroscopic quantities, we consider an arbitrary 
bmction T{S'\S) that satisfies 

T{S'\S)W{S) = T{S\S')W{S'). (6) 

While T{S'\S) could be the transition probability from S 
to S", in which case the above equation is nothing but the 
detailed balance condition, it does not have to be so for 
the derivation of the BH relation. We sum up each side 
of the above equation over such S and S' that n{S) — n 
and n{S') = n' , yielding 

{T{n'\S))^f{n)oj{n) = {T{n\S'))^, f{n')u:{n') (7) 

where 

T{n'\S)^ ^ T{S',S), 
s' 

(„(S')=n') 

and (• • •)„ is the microcanonical average, i.e., 
(Q(5))„^ ^ QiS)^ ^ 1. 

(,i(S) = 7i) {'^(S) = 7^) 

For the derivation of the BH relation, we only have to 
consider the case where S' and S have the same spin 
configuration in common. Then, the procedure for 
generating 5' = (^',F') from S = (*,F) is as follows: 

1. Generate an integer / (1 < ? < uniform-randomly. 

2. If 7; — 0, generate an integer b {I < b < M) uniform- 

randomly, and change F into F' = F U {w} (i.e., 
activate v) with the probability pt{S) where v = 
il,b). 

2'. If 7; = 1, change F into F' = F\{w} (i.e., inactivate 
v) with the probability Py (S) where v is the vertex 
(/, b) such that 'yi^t — 1. 

3. Repeat 1 and 2 (or 2') LM times. 
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The overall transition probability corresponding to one 
repetition of 1 and 2 (or 2') can be written in the form 
of the transition matrix as 

TiS'\S) = l.J2(^5,,,oJ2j^PtbiS)A{r\TU{{lb)}) 

I b 

+^7.1 E K.mbiS)Mr', r\{{ib)})) , 



where some Kronecker's delta's have been written as 
A{A, B) instead of Sa.b for clarity. In order for this 
transition probability to satisfy © with ISJ and Q), 
must satisfy 



1 

M 



where n = n[S) and v is the vertex at which S and S' 
differ from each other. The simplest choice is 

pt{S)= / \( . and p- 5 = A vl/„ . 
L-n f{n) 

(8) 

Note here that we do not have to make p^ smaller than 
1 if we only need to make T{S'\S) satisfy I© and if T 
does not have to be a real transition matrix. Since this 
is the case for obtaining the BH relation, we shall forget 
about the condition pf<l for a time being. With (jSJ, 
T{n + 1\S) and T{n - 1\S) become 



Tin+l\S) = jJ26,,oJ2 



L-n f{n) 



and 



T{n-l\S) = -^E'^7!,iE'^'ri,i>,i^(*«') = 7;(»T--"-kink), 



where rikink is the number of kinks, i.e., the vertices at 
which A(\['„) = 0. Thus we obtain from |7J for n' — n+1, 



L — n 



uj{n) = {n + 1- (rikink) ri+ 1 ) ^^(n + !)■ 

(9) 

We call this relation the BH relation for quantum sys- 
tems. When multiplied U!{n + l)/uj{n) for n = 0, 1, 2, • • •, 
this yields 



Lu{n) = 



LI 



2-n 



m - (rikink) m, 



(10) 



where uj{0) — 2^ has been used. Note that the validity 
of ® and l(Tn|l does not depend on what transition prob- 
ability we use in the simulation although we have used a 
particular form of it in deriving the relation ^ . There- 
fore, we can use any algorithm to estimate the DOS as 
long as it produces the correct microcanonical ensemble. 

In fact, since we have not guaranteed that < 1 
above in JSJ, T{S'\S) mentioned above may not serve as 



the transition probability. To obtain the real transition 
probability, we simply replace (jSJ by 

pt{S) = min(l, Rn(s)) and 
p-(5) = min(l,i?-ig)_JA(vI/,), 



with 



Rr, 



Mu{^,) f{n + 1) 



L-n f{n) 

This is the transition probability that is actually used in 
the simulation presented below. 

THE SECOND PHASE 

We now turn to the second phase in which the state 
S = (^',r) is updated to 5" = («'',r) with T fixed. This 
can be done by either of the loop algorithm or the worm- 
type algorithm. We here adopt the loop type update [isj . 
The weight with which '5 must be sampled is W{'^,T). 
Following the general framework of the loop algorithm 
we decompose the weight at every active vertex w in F as 



E 



a(a)A(*.,G,) 



(11) 



where Gy is a graph and A(\I'„, G„) is the function that 
takes on or 1 depending on the matching between 
and Gy. The readers are referred to Ref.[l9( for the de- 
tails of the graphs and the A-functions. Here we only 
show an example for the compactness of the description. 
In the case of the 5 = 1/2 antiferromagnetic Heisenberg 
model, the decomposition consists of a single term since 
the local Hamiltonian can be expressed as a single delta 
function. 

u{^y)=-^A{^y,Gf) 

where G^ is a horizontal graph that consists of two 
horizontal lines connecting (i,Z) to (j, Z) and + 1) 
to + 1). The function A(^',„,G'^) is defined for 

V = as 



A(*.,G«: 



1 (M^) + Mj) = V'i+i(i) + ^i+iU) = 0) 

(Otherwise) 



Here we have assumed that the state jr/;;) is the si- 
multaneous eigenstate of all the z-components of spins. 
Namely, tpi{i) {= ±1/2) denotes the eigenvalue of Sf of 
the eigenstate \tpi). Therefore, the update is done by 
placing horizontal graphs on all the vertices in F, and 
connect the end points of horizontal lines by vertical lines 
to form loops, then finally flip every loop with probabil- 
ity 1/2. The generalization to the case where the expan- 
sion of the Hamiltonian (|ll|l consists of multiple terms 
is straight-forward. In such a case, we simply choose a 
graph from the ones for which A(^'„, G^) = 1 with the 
probability proportional to a{Gy). 
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RESULTS 

In order to see if the present procedure yields better 
statistics than previous methods, we calculate the spin- 
1/2 antiferromagnetic Heisenberg chain of ten sites. We 
set the initial value of the modification factor to e~^^^ for 
the whole calculation in the paper. (For an appropriate 
initial value of the modification factor, see the discussions 
in the references^, ^3].) 

We set the cut off L to be 500. This limits the acces- 
sible temperature range to T > 0.05 J. In the first stage, 
it takes approximately 6 x 10^ Monte Carlo steps (MCS) 
before the last run with — InA^ < 10^^ terminates. In 
the second stage of simulation, we fix f{n) and let the 
random walker travel for 1 x 10^ MCS. The final estimate 
of the DOS is shown in Fig. By using the DOS in 
Fig. 12 canonical (fixed-temperature) averages of the en- 
ergy, the specific heat, the free energy, and the entropy 
are calculated as functions of temperature. 

In Fig. 13 we show the relative error in the free energy. 
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The top solid line is the one that is obtained based on 
the working estimate of w(n) resulting from the first stage 
that takes about 6 x 10^ sweeps. The bottom dashed line 
represent the one based on the final estimate of Lo{n) ob- 
tained from the second stage using the estimator (|10|1 . 
The middle dotted line also represent the free energy 
based on the second stage data but without using the 
estimator This corresponds to the standard proce- 

dure of the multicanonical simulation 0, 0| • Namely, the 
middle curve is based on the naive estimate of the DOS 
using the histogram itself, i.e., the estimate using the 
relation Lo{n) cx iJ(n)(jj*''"*'(n) where w*'''^' cx {J{n))^^ is 
the trial DOS estimated in the first stage and H{n) is the 
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FIG. 2: The density of states for the series orders for the 
antiferromagnetic spin-1/2 Heisenberg chain on ten-site. The 
cut off of series orders is L = 500. The first stage of simulation 
was performed for 6 x 10^ Monte Carlo sweeps and the second 
stage was calculated for 1 x 10* Monte Carlo sweeps. 



number of visits to the value n counted during the sec- 
ond stage simulation. The difference between the middle 
curve and the bottom one is as much as 1 to 2 orders of 
magnitude, and it clearly shows the utility of the estima- 
tor 11U|1 . The horizontal dashed line indicates the level 
of the largest error of the previous simulation by Troyer, 
Wessel and Alet0 for the same system with the same 
size. The number of Monte Carlo sweeps performed in 
their simulation was the same as the second stage of our 
simulation. It is natural that their result is about the 
same position as the worst (i.e., the highest) point in the 
middle curve, since the last few runs in the WL method, 
where the modification factor is very close to one, are 
almost equivalent to the ordinary multicanonical simula- 
tion with fixed weight. 

In order to see the system-size dependence of the effi- 
ciency, we calculate the standard deviation a of the spe- 
cific heat per site at temperature T = 0.5J for various 
system sizes. Again, we compare three results: (1) the 
one based on a WL method only with the termination 
condition — InAi < 10~^ (denoted as "Wang-Landau" in 
the figure), (2) the one based on the two stage simula- 
tion with the termination condition — In < 10~^ for 
the first stage and with the naive estimation of the DOS 
using the histogram itself in the second stage (denoted 
as "multicanonical"), and (3) the same as the above but 
with the estimator (|10() (denoted as "broad histogram"). 
In (2) and (3), the total number of Monte Carlo sweeps 
are chosen to be the same as that is spent in (1) in order 
to make the comparison fair. We choose L and M so 
that L/M — 50 is fixed. With this choice the accessible 
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FIG. 3: Relative errors of the free energy for the antifer- 
romagnetic spin-1/2 Heisenberg chain on ten-site. The cut 
off of series orders is L = 500. The solid line is the result 
of the Wang-Landau algorithm of the first stage, which took 
6 X 10® MCS. The dotted line represents the result of the mul- 
ticanonical method in the second stage of simulation. The 
short dashed line is the result by using the broad histogram 
relation in the second stage of simulation. The dashed line 
is the result of Troyer et. al. 9] that is drawn at the worst 
value. 
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temperature range is T > 0.05 J. 40 independent sets are 
performed for M = 8, 12, 16, 20 whereas 12 for M = 24. 
The results are shown in Fig. ^ The system size de- 
pendence is small. The results of the WL algorithm and 
the naive estimate of the DOS based on the two-stage 
simulation are close to each other again. The proposed 
method outperforms the other method by more than an 
order of magnitude in all the cases studied here. 

SUMMARY 

Generalizing the histogram methodsll H E3 

based 

on the graph representation, we have presented a Monte 
Carlo method for quantum lattice models that consists 
of two stages. For the method's performance, the gener- 
alized broad histogram relation 0, proved for quantum 
systems, is crucial. In the first stage, the WL method is 
employed as in 9] for obtaining the working estimate of 
the DOS with a short CPU time. In the second stage, 
microcanonical averages are computed and the final esti- 
mate of the DOS is obtained from them through the use 
of the generalized BH relation. In both the stages, the 
update of the states are performed by loop algorithm. We 
have demonstrated that the proposed method yields very 
accurate estimates of the DOS in the case of the = 1/2 
antiferromagnetic Heisenberg model, which we believe to 
be a typical case that represents many other cases. In 
Ref. [3], an alternative usage of the extended ensemble 
method was proposed in which one takes the series ex- 




o Wang-Landau 
□ multicanonical 
o broad histogram 




FIG. 4: The system-size dependences of the standard de- 
viation (7 for the specific heat per one site at T = 0.5J 
for the spin-1/2 antiferromagnetic Heisenberg chain. Three 
methods were compared: the generalized WL algorithm (cir- 
cles), the combination of the multicanonical method and 
the generalized WL algorithm (squares), and the proposed 
method (diamonds). The methods spent the same MCS 
for each M. We took the M and the L as (M, L) = 
(8, 400), (12, 600), (16, 800), (20, 1000). The accessible tem- 
perature range is T > 0.05J. 40 sets of simulation are per- 
formed for M = 8, 12, 16, 20 and 12 sets for M = 24. 



pansion not in /3 but in other coupling constants. This 
type of extension is also available in the present method. 
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